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Abstract 

We present a novel approach for improving particle filters for multi- 
target tracking. The suggested approach is based on drift homotopy 
for stochastic differential equations. Drift homotopy is used to design 
a Markov Chain Monte Carlo step which is appended to the particle 
filter and aims to bring the particle filter samples closer to the observa- 
tions. Also, we present a simple Metropolis Monte Carlo algorithm for 
tackling the target-observation association problem. We have used the 
proposed approach on the problem of multi-target tracking for both 
linear and nonlinear observation models. The numerical results show 
that the suggested approach can improve significantly the performance 
of a particle filter. 



Introduction 



Multi-target tracking is a central and difficult problem arising in many sci- 
entific and engineering applications including radar and signal processing, 
air traffic control and GPS navigation [llj. The tracking problem consists 
of computing the best estimate of the targets' trajectories based on noisy 
measurements (observations) . 

Several strategies have been developed for addressing the multi-target 
tracking problem, see e.g. [Il[6l[5l[l4l|8l[l0l[IIl[l2l[l3l[l9]. In this paper we 
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focus on particle filter techniques OdJj. The popularity of the particle filter 
method has increased due to its flexibility to handle cases where the dynamic 
and observation models are non- linear and/or non-Gaussian. The particle 
filter approach is an importance sampling method which approximates the 
target distribution by a discrete set of weighted samples (particles). The 
weights of the samples are updated when observations become available in 
order to incorporate information from the observations. 

Despite the particle filter's flexibility, it is often found in practice that 
most samples will have a negligible weight with respect to the observation, 
in other words their corresponding contribution to the target distribution 
will be negligible. Therefore, one may resample the weights to create more 
copies of the samples with significant weights |8j. However, even with the 
resampling step, the particle filter might still need a lot of samples in order 
to approximate accurately the target distribution. Typically, a few samples 
dominate the weight distribution, while the rest of the samples are in sta- 
tistically insignificant regions. Thus, some authors (see e.g. [3 [20]) have 
suggested the use of an extra step, after the resampling step, which can help 
move more samples in statistically significant regions. 

The extra step for the particle filter is a problem of conditional path 
sampling for stochastic differential equations (SDEs). In [16], a new ap- 
proach to conditional path sampling was presented. In that paper, it was 
also shown how the algorithm can be used to perform the extra step of a 
particle filter. In the current work, we have applied the conditional path 
sampling algorithm from |16j to perform the extra step of a particle filter 
for the problem of multi-target tracking. 

The suggested approach is based on drift homotopy for the SDE system 
which describes the dynamics of the targets. The dynamics of an SDE are 
governed by a deterministic term, called the drift, and a stochastic term, 
called the noise. While unconditional path sampling is straightforward for 
SDEs, albeit expensive for high dimensional systems, conditional path sam- 
pling can be difficult even for low dimensional systems. On the other hand, 
it can be easier to find conditional paths for an SDE with a modified drift 
which is usually simpler than the drift of the original equation. Of course, 
these simplified paths may have a very low probability of being paths of 
the original SDE. Drift homotopy proceeds by considering a sequence of 
SDEs with drifts which interpolate between the original and modified drifts. 
Then one samples (through MCMC sampling) paths from each SDE in the 
sequence using the best sample from the previous SDE in the sequence as 
the initial condition for the MCMC sampling. This allows one to gradu- 
ally morph a path of the modified SDE (which maybe easier to satisfy the 
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conditions) to a conditional path of the original SDE. 

In addition to the extra step for the particle filter, we have designed 
and implemented a simple Metropolis Monte Carlo algorithm for the target- 
observation association problem. This algorithm effects a probabilistic search 
of the space of possible associations to find the best target-observation as- 
sociation. Of course, one can use more sophisticated association algorithms 
(see jl4] and references therein) but the Monte Carlo algorithm performed 
very well in the numerical experiments. 

This paper is organized as follows. Sections 11.11 and 11.21 provide a brief 
presentation of particle filters for single and multiple targets (more details 
can be found in [8l [H [U [H]), which will serve to highlight the versatil- 
ity and drawbacks of this popular filtering method. Sections 11.31 and 11.41 
demonstrate how one can use an extra step to improve the performance of 
particle filters for single and multiple targets. In particular, we discuss how 
drift homotopy can be used to effect the extra step. Section [2] describes 
the Monte Carlo sampling algorithm for computing the best association be- 
tween observations and targets for the case of multiple targets. Section [3] 
presents numerical results for multi-target tracking for the cases of linear 
and nonlinear observation models. Finally, Section [J] contains a discussion 
of the results as well as directions for future work. 

1 Particle filtering 

Particle filters are a special case of sequential importance sampling methods. 
In Sections 11.11 and 11.21 we discuss the generic particle filter for a single and 
multiple targets respectively. In Sections 1 1.31 and II. 41 we discuss the addition 
of an extra step to the generic particle filter for the cases of a single and 
multiple targets respectively. 

1.1 Generic particle filter for a single target 

Suppose that we are given an SDE system and that we also have access 
to noisy observations Zt^ , ■ ■ ■ , Ztj^ of the state of the system at specified 
instants Ti, . . . , Tk- The observations are functions of the state of the sys- 
tem, say given by Zt^. = G{XT^,^k), where ^k,k = 1,...,K are mutually 
independent random variables. For simplicity, let us assume that the distri- 
bution of the observations admits a density g(XT,., Zt,.), i.e., p{ZTf.\XT,.) oc 
c,(Xt,,ZtJ. 

The filtering problem consists of computing estimates of the conditional 
expectation E[f{XT^.)\{Z'rj}j^i], i.e., the conditional expectation of the 
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state of the system given the (noisy) observations. Equivalently, we are 
looking to compute the conditional density of the state of the system given 
the observations p{XT^.\{ZTj}j=i)- There are several ways to compute this 
conditional density and the associated conditional expectation but for prac- 
tical applications they are rather expensive. 

Particle filters fall in the category of importance sampling methods. Be- 
cause computing averages with respect to the conditional density involves 
the sampling of the conditional density which can be difficult, importance 
sampling methods proceed by sampling a reference density q{XT^.\{ZTj}j=i) 
which can be easily sampled and then compute the weighted sample mean 



or the related estimate 



where N has been replaced by the approximation 



f^q{Xl\{Zr,}U) 



Particle filtering is a recursive implementation of the importance sampling 
approach. It is based on the recursion 

p{XtJ{Zt^}';=^) ^ giXT„ZTMXTj{ZT,}';zl), (2) 
where p(XTj{ZT,},t!) = |p(XtJXt,_ Jp(Xt,_ J{Zt, )dXT,_, . (3) 

If we set 

q{XTj{ZT^}U)=PiXn\{ZT,}';Zl), 
then from ([2]) we get 



p{Xt,\{Zt^}';=,] 

q{XT,\{ZTj^^,] 



(xg{XT„ZT,). 
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The approximation in expression ([T]) becomes 



E[f{XT.)\{ZTAU - I' \ (4) 



From Q we see that if we can construct samples from the predictive distri- 
hiitinn TtfYrr. \{:^Tj}'^zl) then we can define the (normalized) weights = 



g{XJ^ ,Zt.) 

7^- — - — -, use them to weigh the samples and the weighted samples 

will be distributed according to the posterior distribution pC-^^t^ |{^Tj }j=i)- 
In many applications, most samples will have a negligible weight with 
respect to the observation, so carrying them along does not contribute sig- 
nificantly to the conditional expectation estimate (this is the problem of 
degeneracy [9]). To create larger diversity one can resample the weights to 
create more copies of the samples with significant weights. The particle filter 
with resampling is summarized in the following algorithm due to Gordon et 

ai. [g. 

Particle filter for a single target 

1. Begin with unweighted samples X^^ ^ from p(X2^j._^ [{Zt^. 

2. Prediction: Generate samples irom p{XTf.\XTi._j^)- 

3. Update: Evaluate the weights 

T.n=l 9{X'^,^ZtS 

4. Resampling: Generate N independent uniform random variables 
{6'"}^^i in (0, 1). Torn = 1, . . . , let = where 

1=1 1=1 
where j can range from 1 to N. 

5. Set k = k + 1 and proceed to Step 1. 

The particle filter algorithm is easy to implement and adapt for dif- 
ferent problems since the only part of the algorithm that depends on the 
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specific dynamics of the problem is the prediction step. This has led to the 
particle filter algorithm's increased popularity [5]. However, even with the 
resampling step, the particle filter can still need a lot of samples in order to 
describe accurately the conditional density p(-'^Tfe K-^T, }j=i)- Snyder et al. 
[15| have shown how the particle filter can fail in simple high dimensional 
problems because one sample dominates the weight distribution. The rest 
of the samples are not in statistically significant regions. Even worse, as we 
will show in the numerical results section, there are simple examples where 
not even one sample is in a statistically significant region. In the next sub- 
section we present how drift homotopy can be used to push samples closer 
to statistically significant regions. 

1.2 Generic particle filter for multiple targets 

Suppose that we have A = 1, . . . , A targets. Also, for notational simplicity, 
assume that the Ath target comes from the Ath observation. Even when this 
is not the case, we can relabel the observations to satisfy this assumption. 
We will discuss in Section [2] how targets can be associated to observations. 
The targets are assumed to evolve independently so that the observation 
weight of a sample of the vector of targets is the product of the individual 
observation weights of the targets [T^. The same is true for the transition 
density of the vector of targets between observations. We denote the vector 
of targets at observation Tk by 

= (^i,Tfe, • • • ,^A,rJ 
and the observation vector at by 

Also, we can have different observation weight densities g\, A = 1, . . . , A for 
different targets. However, in the numerical examples we have chosen the 
same observation weight density for all targets. 

Following |14] we can write the particle filter for the case of multiple 
targets as 

Particle filter for multiple targets 

1. Begin with N unweighted samples from p{XTf._^\{ZTj}^zl) = 



6 



2. Prediction: Generate samples from 

A 

p{XtJXt,_,) = llp{Xx,n\Xx,n.,). 

A=l 

3. Update: Evaluate the weights 

^ Ux=i9x{X'x,n^Zx,n) 

Yln=lU\=l9x{X'l^Tk^Z>;n) 

4. Resampling: Generate N independent uniform random variables 
{e''}n=i in (0, 1). For n = 1, . . . , let = where 

1=1 1=1 
where j can range from 1 to N. 

5. Set k = k + 1 and proceed to Step 1. 

1.3 Particle filter with MCMC step for a single target 

Several authors (see e.g. [3 [20]) have suggested the use of a MCMC step 
after the resampling step (Step 4) in order to move samples away from 
statistically insignificant regions. There are many possible ways to append 
an MCMC step after the resampling step in order to achieve that objective. 
The important point is that the MCMC step must preserve the conditional 
density p{XT^.\{ZTj}j=l)■ In the current section we show that the MCMC 
step constitutes a case of conditional path sampling. 

We begin by noting that one can use the resampling step (Step 4) in the 
particle filter algorithm to create more copies not only of the good samples 
according to the observation, but also of the values (initial conditions) of the 
samples at the previous observation. These values are the ones who have 
evolved into good samples for the current observation (see more details in 
[20]). The motivation behind producing more copies of the pairs of initial 
and final conditions is to use the good initial conditions as starting points 
to produce statistically more significant samples according to the current 
observation. This process can be accomplished in two steps. First, Step 4 
of the particle filter algorithm is replaced by 



7 



Resampling: Generate independent uniform random variables {0""}^ 
in (0, 1). Forn = 1, . . . , TV let (X^^_^,X^J = (X^^_^,X^Jwliere 

1=1 1=1 

Also, through Bayes' rule [20] one can show that the posterior density 
|{^Tj }j=i) is preserved if one samples from the density 

g{XT^,ZT^)p{XT^\XT^_^) 

where ^t^.i are given by the modified resampling step. This is a problem 
of conditional path sampling for (continuous-time or discrete) stochastic 
systems. The important issue is to perform the necessary sampling efficiently 

We propose to do that here using drift homotopy. In particular, suppose 
that we are given a system of stochastic differential equations (SDEs) 

dXt = a{Xt)dt + a{Xt)dBt, (5) 

Also consider an SDE system with modified drift 

dYt = b{Yt)dt + a{Yt)dBt, (6) 

where b{Yt) is suitably chosen to facilitate the conditional path sampling 
problem. 

Consider a collection of L + 1 modified SDE systems 

dy/ = (1 - ei)b{Y})dt + €ia{Y})dt + cj{Yl)dBt, 

where G [0, 1], / = 0, . . . , L, with e/ < e^+i, eo = and = Instead of 
sampling directly from the density g{XT,. , Zt^. )p{Xt^, \Xt,._^) , one can sample 
from the density giYj,^^ ZT^)p{Y^J\XTf,_^) and gradually morph the sample 
into a sample of g{XT^^ , Zt^^ )p(^Tfe l-^^T^.i ) • 

Drift homotopy algorithm: 

• (/ = 0) Begin with a sample from the modified SDE Q. 

• Sample through MCMC the density g{Y^ ,ZTf?jp{Y:^ \^n_i)- 
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• For / = 1, L take the last sample from the (/ — l)st SDE and use it 
as in initial condition for MCMC sampling of the density 

at the Ith level. 

• Keep the last sample at the Lth level. 

The drift homotopy algorithm is similar to Simulated Annealing (SA) used 
in equilibrium statistical mechanics [DJ. However, instead of modifying a 
temperature as in SA, here we modify the drift of the system. 

We are now in a position to present the particle filter with MCMC step 
algorithm 

Particle filter with MCMC step for a single target 

1. Begin with N unweighted samples X^^ ^ from 

2. Prediction: Generate samples from p{Xt^.\Xt^._j^)- 

3. Update: Evaluate the weights 

4. Resampling: Generate N independent uniform random variables 
{^^=1 in (0,1). For n = 1,...,A let {X^^_^,X^J = {X'l_^,X'l) 
where 

1=1 1=1 

where j can range from 1 to N. 

5. MCMC step: For n = 1,...,A choose a modified drift (possibly 
different for each n). Construct a path for the SDE with the modified 
drift starting from XJj^^ ^ . Construct through drift homotopy a Markov 
chain for with stationary distribution 

6. SetX" =y". 
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7. Set k = k + 1 and proceed to Step 1. 



Since the samples X^^ = are constructed by starting from different 
sample paths, they are independent. Also, note that the samples Xj,^ are 
unweighted. However, we can still measure how well these samples approx- 
imate the posterior density by comparing the effective sample sizes of the 
particle filter with and without the MCMC step. For a collection of A'^ 
samples the effective sample size ess{Tk) is defined by 

ess{Tk) = 



where 



r - ^ 



N N 

-Y^[g[XJ},^,ZT,)-W,)^ and W, = -Y.9{Xn^ZT,). 

n=l n=l 



The effective sample size can be interpreted as that the N weighted samples 
are worth of ess{Tk) = i.i.d. samples drawn from the target density, 

which in our case is the posterior density. By definition, ess{Tk) < N. If the 
samples have uniform weights, then ess{Tk) = N. On the other hand, if all 
samples but one have zero weights, then ess{Tk) = 1. 



1.4 Particle filter with MCMC step for multiple targets 

We discuss now the case of multiple, say A, targets. Instead of the observa- 
tions for a single target now we have a collection of observations for all the 
targets {Zt,%^ = {{Z\^.. . . , 

The particle filter with MCMC step for the case of multiple targets is 

Particle filter with MCMC step for multiple targets 

1. Begin with A unweighted samples XJy^ ^ from p(Aj'^_j }^^|) = 
nA=lP(^A,T,_J{^A,T,},t!). 

2. Prediction: Generate A samples Aj? from 

A 

p(AtJAt,_J = nP(^A,TjA,,T,_J. 
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3. Update: Evaluate the weights 



En=i 111=1 9\ {x'Xn ' ^^."^fc ) 

4. Resampling: Generate N independent uniform random variables 
{n;r=i in (0,1). For n=\,...,N let = {X'l_^,X'l) 
where 

1=1 1=1 

where j can range from 1 to N. 

5. MCMC step: For n = 1, . . . ,N and A = 1, . . . , A choose a modi- 
fied drift (possibly different for each n and each A). Construct a path 
for the SDE with the modified drift starting from X^j,^_^. Construct 
through drift homotopy a Markov chain for with stationary dis- 
tribution 



n5A(iT,^A,TjPA(lTI^AV._,)- 



A=l 

6. SetX^^=Kp^. 

7. Set k = k + 1 and proceed to Step 1. 

For a collection of N samples the effective sample size essA(Tfc) for A 
targets is 

essA{Tk] 



where 



CAh 



N A 



WAk\ iv^^J-J-- 

^^''^ \ n=l A=l 

^ AT A 

and WA,k = E n 5a(^aV,, ^a,tJ. 



n=l A=l 
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2 Monte Carlo target-observation association al- 
gorithm 



As in any other sequential Monte Carlo algorithm for multi-target tracking 
(see e.g. [H] and references therein) we need to associate, for each sam- 
ple, the evolving targets to the observations. The association problem is 
sensitive to the tracking accuracy of the algorithm. If we cannot follow 
accurately each target and two or more targets are close, then the associ- 
ation algorithm can assign the wrong observations to the targets. After a 
few observation steps this can lead to the inability to follow a target's true 
track anymore. There are different ways to perform the target-observation 
association. The ones that we are aware of are based on various types of as- 
signment algorithms first developed in the context of computer science |3l [2] . 
Here we have decided on using a different algorithm to perform the target- 
observation association. In particular, we have designed a simple Metropolis 
Monte Carlo algorithm which effects a probabilistic search of the space of 
possible associations to find the best target-observation association. 

Suppose that at observation time t we have A^ observations which include 
surviving and possible newborn targets. Also, suppose that for each target 
we have at each step J observation values for different quantities depending 
on the target's state. For simplicity, we assume that we observe the same 
quantities for all targets. Usually, one observes position related quantities, 
like the x, y position or the bearing and range of the target. The observation 
model is given by 

Zx,j,t = gji^x,t) + Caj, for A = 1, . . . , At and j = 1,...,J, (7) 

where for A = 1, . . . , A^ and j = 1, . . . , J, are i.i.d. random variables. 
For simplicity let us assume that the S^xj are ~ N{0, c7|). Note that the ob- 
servation model ([7]) assumes that we know that the Ath observation comes 
from the Ath target. However, in reality, we do not have such information 
and we need to make an association between observations and targets. De- 
fine an association map A given by A(X) = mx where X,mx = 1,... ,At. 
The association map assigns to the Ath observation the m^^th target. For 
At observations there are A^! different observation-target association maps. 
For a specific association of each observation to some target, the likelihood 
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function for the collection of observations for the nth sample is given by 



oc ji n 

A=li=l 



(8) 



where the proportionality is up to an (immaterial for our purposes) normal- 
ization constant. If we define 



Hobsi^U^ ■ ■ ■ ,^At,tiZl,l,t, ■ ■ ■ ,Zi,j,t, ■ ■ ■ ,ZAt,J,t) 

(^A,,..-^.K(A).t))' 

we can rewrite p (omitting all the arguments of p and for simplicity) 
as 

pcxexp[-Hi,] (9) 

The superscript A is to denote the dependence of the value of H^f^^ on 
the specific association map A. The best association map is the one which 
maximizes p. By definition, > and thus the best association is the one 
which minimizes -ff^^- We can find the best association map by standard 
Metropolis Monte Carlo sampling on the space of association maps, using 
the density p (see e.g. [9j about Metropolis Monte Carlo sampling). 

As in every Metropolis Monte Carlo sampling algorithm there is arbi- 
trariness in the new configuration (association map in our case) proposal. 
We have tried two different proposal schemes which performed equally well, 
at least for cases up to 4 targets that we examined in our numerical ex- 
periments. The first proposal scheme constructs a new observation-target 
association map from scratch. This means that for each observation in the 
observation vector we choose a new target to associate it with. Note that if 
there is only one target there is no need for sampling since there is only one 
possible observation-target association. The second proposal scheme starts 
with an initial association map. Then it chooses randomly a pair of obser- 
vations and their associated targets and exchanges the associated targets. 
This allows one to build a good association map incrementally. We expect 
that the second proposal scheme will be advantageous in the case when 
the number of observations at the observation instant t is large. Recall 
that the number of possible association maps grows as A^!, and it becomes 
prohibitively large for an exhaustive search even for moderate values of Aj. 
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The numerical results we present here are for the first proposal scheme. 
The Metropolis sampling algorithm was run with 10000 Metropolis ac- 
cept/reject steps and we kept the last accepted association map. 



3 Numerical results 

We present numerical results for multi-target tracking using the particle 
filter with an MCMC step performed by drift homotopy and hybrid Monte 
Carlo. We have synthesized tracks of targets moving on the xy plane using 
a 2D near constant velocity model [1]. At each time t we have a total of At 
targets and the evolution of the Ath target (A = 1, . . . , A^) is given by 



XA,t = AxA,t-l + BvA,t 

= [x\,t,xx,t,y\,t,y\,tV , 



(10) 



where (xA,t,XA,t) and [y\^f,y\,t) are the xy position and velocity of the Ath 
target at time t. The matrices A and B are given by 



1 


T 











1 














1 


T 











1 



and B 



" rV2 





T 













T 



(11) 



where T is the time between observations. For the experiments we have set 
T = 1, i.e., noisy observations of the model are obtained at every step of the 
model (jlOp . The model noise 'v^ t is a collection of independent Gaussian 
random variables with covariance matrix I]„ defined as 







(12) 



In the experiments we have o"^ ~ ^"y ~ 1- Also, we have considered two 
possible cases for the observation model, one linear and one nonlinear. Due 
to the different possible combinations of targets to observations we use a 
different index m to denote the obsevations. Since we do not assume any 
clutter we have m = 1, . . . , A^. If the mth observation z^.t at time t comes 
from the Ath target we have 



yu 



+ w 



(13) 
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for the linear observation model and 



arctan 



+ yit) 



y\,t \ 

2 \l/2 



+ w. 



(14) 



for the nonlinear observation model. As is usual in the literature, the non- 
linear observation model consists of the bearing 6 and range r of a target. 
The observation noise Wm t is white and Gaussian with covariance matrix 



<hs,x 



a 



obs,y 



for the linear observation model and 



ag U 




(15) 



(16) 



for the nonlinear observation model. For the numerical experiments with 



the linear observation model we chose 0"^^^ ^ 



obs,y 



1. For the numerical 



experiments with the nonlinear observation model we chose = 10~^ and 
(T^ = 1. These values make our example comparable in difficulty to examples 
appearing in the literature (see e.g. [HI ITSl 119] ). 

The synthesized target tracks were created by specifying a certain sce- 
nario, to be detailed below, of surviving, newborn and disappearing targets. 
According to this scenario we evolved the appropriate number of targets 
according to (fTOj) and recorded the state of each target at each step. For 
the surviving targets we created an observation by using the state of the 
target in the observation model. Thus, for the linear observation model, 
the observations were created directly in xy space by perturbing the xy 
position of the target by (jl3p . For the nonlinear observation model, the ob- 
servations were created in bearing and range space 9,r by using ()14p . The 
perturbed bearing and range were transformed to xy space by the trans- 
formation X = rcos6, y = rsinO to create a position for the target in xy 
space. 

The newborn targets for the linear model were created in xy space di- 
rectly by sampling uniformly in [—100, 100]. Afterwards, the observations of 
the newborn targets were constructed by perturbing the x, y positions using 
(jl3p . The newborn targets for the nonlinear model were created in xy space 
by sampling uniformly in [—100, 100]. Afterwards, we transformed the x,y 
positions to the bearing and range space 9, r and perturbed the bearing and 
range according to ()14p . The perturbed bearing and range were again trans- 
formed back to xy space to create the position of the newborn target. Note 
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that both observation models do not involve the velocities. The newborn 
target velocities were sampled uniformly in [—1, 1]. 

The number of targets at each observation instant is: Aq = 2, Ai = 2, 
A2 = 1, A3 = 2, A4 = 3, Ag = Ag = . . . = A200 = 4. So, for the majority 
of the steps we have 4 targets which makes the problem of tracking rather 
difficult. 



3.1 Drift homotopy 

The dynamics of the targets for the modified drift systems are given by 

^{,t = ^At-i + + BvA,t, 

where Z\^\^fi Zz,\,t ^2,\,ti ZA,\,t are the xy positions and velocities respec- 
tively for the Ath target at time t. 



The matrix c' is given by 

= (1 - ei) 



rj-i'2 

fJ'X — 



where € [0, 1], I = 0, . . . , L, with e/ < e/+i, eo = and cl = 1. 

For the nth sample, the density we have to sample for the linear model 

is 

At. 



oc 



JJexp 



A=l 



iZi,x,k - <A,fe-i - ^lx,k-iT - (1 - eOi^l^^ - ^^'^ 



iZ3,x,k - ^3"A,fe-i - ^lx,k-iT - (1 - eOf^l,^ - ^v-,^. 



+ 



obs,y 



+ 



+ 



-'y,\,k) 

2al 



where 



n _ N Tlin'=li^l,\ ,k-l + ^2,A,fe-l^) ^l,X,k-l _ 2^2,A,fc-l 
l^x,X — 



r2/2 
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and 



1 v^A^ , 

_ _, , , ■4,A,fc-l 

^y,>' ~ T2 /9 



3.1.1 Hybrid Monte Carlo 

We chose to use Hybrid Monte Carlo to perform the samphng (any other 
MCMC method can be used) [9j. We present briefly the hybrid Monte Carlo 
(HMC) formulation that we have used to sample the conditional density for 
each target. We will formulate HMC for the Ith level of the drift homotopy 
process. 

Define the potential fc) by 

^ ^ (^l,A,fe - ^r,A,fc-l ~ ^2,A,fc-l^ - (1 - ^OMx.A^ ~ ^^x,A,fc)^ 



E 



A=l oos,x 

^ iZ3Xk - 4,A,fe-i - ^lx,k^iT - (1 - q)/^^,a¥ - ^^lx,k)' 

Define the vectors v^,^ = [v^^^^, <,At,,J^ and v^,^ = [v^^^f^, v^^^^J^^ 
where T is the transpose (not to be confuse with the interval between 
observations T used before). Consider v'^j^^Vy^ as the position variables 
of a Hamiltonian system. We define the 2 At j, -dimensional position vector 
q = [qi, 52]"^ with qi = v"^ and q2 = f.. To each of the position variables 
we associate a momentum variable and we write the Hamiltonian 

H,^{q,p)=V,^{q) + ^, 

where p = [pi,P2]'^ is the momentum vector. Thus, the momenta variables 
are Gaussian distributed random variables with mean zero and variance 1. 
The equations of motion for this Hamiltonian system are given by Hamilton's 
equations 

dqi dH^ dpi dH^ • 1 o 

— = ^ and — = — -— ^ for z = 1, 2. 
ar opi dr oqi 

HMC proceeds by assigning initial conditions to the momenta variables 

T 

(through sampling from exp( — 2^)), evolving the Hamiltonian system in 
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Figure 1: Linear observation model. The solid lines denote the true target 
tracks, the crosses denote the observations and the dots the conditional 
expectation estimates from the MCMC particle filter. We have plotted the 
conditional expectation estimates every 5 observations to avoid cluttering 
in the figure. 

fictitious time r for a given number of steps of size 6t and then using the 
solution of the system to perform a Metropolis accept/reject step (more de- 
tails in [^). After the Metropolis step, the momenta values are discarded. 
The most popular method for solving the Hamiltonian system, which is the 
one we also used, is the Verlet leapfrog scheme. In our numerical imple- 
mentation, we did not attempt to optimize the performance of the HMC 
algorithm. For the sampling at each level of the drift homotopy process we 
used 10 Metropolis accept/reject steps and 1 HMC step of size 6t = 10~^ 
to construct a trial path. A detailed study of the drift homotopy/HMC 
algorithm for conditional path sampling problems outside of the context of 
particle filtering will be presented in a future publication. 

For the nonlinear observation model we can use the same procedure 
as in the linear observation model to define a Hamiltonian system and its 
associated equations. We omit the details. 

3.2 Linear observations 

We start the presentation of our numerical experiments with results for the 
linear observation model ([13]) . Figures [1] and [2] show the evolution in the 
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X 

Figure 2: Linear observation model. Detail of Figure [TJ 

xy space of the true targets, the observations as well as the estimates of 
the MCMC particle filter. It is obvious from the figures that the MCMC 
particle filter follows accurately the targets and there is no ambiguity in the 
identification of the target tracks. 

The performance of the MCMC particle filter with 100 samples is com- 
pared to the performance of the generic particle filter with 120 samples in 
Figure [3] by monitoring the evolution in time of the RMS error per target. 
The RMS error per target (RMSE) is defined with reference to the true 
target tracks by the formula 



RMSE{t) = 




where || • || is the norm of the position and velocity vector. Note that the 
state vector norm involves both positions and velocities even though the 
observations use information only from the positions of a target, x^^t is the 
true state vector for target k. £'[xfc j |Zi, . . . ,Zt] is the conditional expecta- 
tion estimate calculated with the MCMC or generic particle filter depending 
on whose filter's performance we want to calculate. 

The MCMC particle filter has a computational overhead of the order of 
a few percent compared to the generic particle filter. We have thus used 
the generic particle filter with more samples than the MCMC particle filter. 
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Figure 3: Linear observation model. Comparison of RMS error per target 
for the MCMC particle filter and the generic particle filter. 

This additional number of samples more than accounts for the computational 
overhead of the MCMC particle filter. As can be seen in Figure[3]the generic 
particle filter's accuracy deteriorates quickly. On the other hand, the MCMC 
particle filter maintains an 0(1) RMS error per target for the entire tracking 
interval. The average value of the RMS error over the entire time interval of 
tracking is about 2.5 with standard deviation of about 0.5. For the generic 
particle filter, the average of the RMS error over the time interval of tracking 
is about 800 with standard deviation of about 760. 

Figure H] compares the effective sample size for the generic particle filter 
and the MCMC particle filter. Because the number of samples is different 
for the two filters we have plotted the effective sample size as a percentage 
of the number of samples. We have to note that, after about 50 steps, 
the generic particle filter started producing observation weights (before the 
normalization) which were numerically zero. This makes the normalization 
impossible. In order to allow the generic particle filter to continue we chose 
at random one of the samples, since all of them are equally bad, and assigned 
all the weight to this sample. We did that for all the steps for which the 
observation weights were zero before the normalization. As a result, the 
effective sample size for the generic particle filter drops down to 1 sample 
after about 50 steps. Once the generic particle filter deviates from the true 
target tracks there is no mechanism to correct it. Also, we tried assigning 
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Figure 4: Linear observation model. Comparison of effective sample size for 
the MCMC particle filter and the generic particle filter. 

equal weights to all the samples when the observation weights dropped to 
zero. This did not improve the generic particle's performance either. On 
the other hand, the MCMC particle filter maintains an effective sample size 
which is about 25% of the number of samples. 

3.3 Nonlinear observations 

We continue with results for the nonlinear observation model (jl4p . Figures [5] 
and [6] show the evolution in the xy space of the true targets, the observations 
as well as the estimates of the MCMC particle filter. Again, as in the case 
of the linear observation model, the MCMC particle filter follows accurately 
the targets and there is no ambiguity in the identification of the target 
tracks. 

The case of the nonlinear observation model is much more difficult than 
the case of the linear observation model. The reason is that for the nonlinear 
observation model, the observation errors, though constant in bearing and 
range space, they become position dependent in xy space. In particular, 
when X and/or y are large, the observation errors can become rather large. 
This is easy to see by Taylor expanding the nonlinear transformation from 
bearing and range space to xy space around the true target values. Suppose 
that the true target bearing and range are ^Oj^o and its xy space position 
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Figure 5: Noninear observation model. The solid lines denote the true target 
tracks, the crosses denote the observations and the dots the conditional 
expectation estimates from the MCMC particle filter. We have plotted the 
conditional expectation estimates every 5 observations to avoid cluttering 
in the figure. 
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Figure 6: Nonlinear observation model. Detail of Figure [5l 
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is xq = rQCOs9Q,yQ = rosin^o- Also, assume that the observation error in 
bearing and range space is, respectively, 59 and 5r. The xy position of a 
target that is perturbed by 69 and Sr in bearing and range space is (to first 
order) 

X = xq — yo59 — 5r cos 
y = yo + xq59 — 6r sin ^o- 

Thus, the perturbation in xy space can be significant even if 59 and 5r are 
small. In our example we have = 10~^. So, when the true target x and 
y values become of the order of 10^ as happens for some of the targets, the 
observation value in bearing and range space can be quite misleading as far 
as the xy space position of the target is concerned. As a result, even if one 
does a good job in following the observation in bearing and range space, the 
conditional expectation estimate of the xy space position can be inaccurate. 

With this in mind, we have used 200 samples for the MCMC particle 
filter and 220 samples for the generic particle filter. Again, the extra samples 
used for the generic particle filter more than account for the computational 
overhead of the MCMC particle filter. The performance of the MCMC 
particle filter is compared to the performance of the generic particle filter in 
Figure [7] by monitoring the evolution in time of the RMS error per target. 
The generic particle filter's accuracy again deteriorates rather quickly. The 
error for the MCMC particle filter is larger than in the linear observation 
model but never exceeds about 80 even after 200 steps when the targets have 
reached large values of x and/or y. The average value of the RMS error over 
the entire time interval of tracking is about 22 with standard deviation of 
about 21. For the generic particle filter, the average of the RMS error over 
the time interval of tracking is about 760 with standard deviation of about 
770. 

Figure [8] compares the effective sample size for the generic particle filter 
and the MCMC particle filter. After about 60 steps, the generic particle fil- 
ter, started producing observation weights (before the normalization) which 
were numerically zero. This makes the normalization impossible. In order 
to allow the generic particle filter to continue we chose at random one of the 
samples, since all of them are equally bad, and assigned all the weight to 
this sample. We did that for all the steps for which the observation weights 
were zero before the normalization. As a result, the effective sample size 
for the generic particle filter drops down to 1 sample after about 60 steps. 
Once the generic particle filter deviates from the true target tracks there is 
no mechanism to correct it. Also, we tried assigning equal weights to all 
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Figure 7: Nonlinear observation model. Comparison RMS error per target 
for the MCMC particle filter and the generic particle filter. 




Figure 8: Nonlinear observation model. Comparison of effective sample size 
for the MCMC particle filter and the generic particle filter. 
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the samples when the observation weights dropped to zero. This did not 
improve the generic particle's performance either. On the other hand, the 
MCMC particle filter maintains an effective sample size which is about 25% 
of the number of samples. 

4 Discussion 

We have presented an algorithm for multi-target tracking which is based on 
drift homotopy for stochastic differential equations. The algorithm builds on 
the existing particle filter methodology for multi-target tracking by append- 
ing an MCMC step after the particle filter resampling step. The purpose of 
the addition of the MCMC step is to bring the samples closer to the obser- 
vation. Even though the addition of an MCMC step for a particle filter has 
been proposed and used before [7j , to the best of our knowledge, the use of 
drift homotopy to effect the MCMC step is novel (see also [16]). 

We have tested the performance of the algorithm on the problem of 
tracking multiple targets evolving under the near constant velocity model [Ij . 
We have examined two cases of observation models: i) a linear observations 
model involving the positions of the targets and ii) a nonlinear observation 
model involving the bearing and range of the targets. For both cases the 
proposed MCMC particle filter exhibited a significantly better performance 
than the generic particle filter. Since the MCMC particle filter requires 
more computations than the generic particle filter it is bound to be more 
expensive. However, the computational overhead of the MCMC particle 
filter is rather small, of the order of a few extra samples worth for the 
generic particle filter. 

We plan to perform a detailed study of the proposed algorithm in more 
realistic cases involving clutter, spawning and merging of targets. Also, we 
want to study the behavior of the algorithm for cases with random birth and 
death events as well as for larger number of targets. Finally, the algorithm 
can be coupled to any target detection algorithm (see e.g. [14]) to perform 
joint detection and tracking. 
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